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Abstract 

A kinetic Monte Carlo (KMC) method is used to study the struc- 
tural properties and dynamics of a supercooled binary Lennard-Jones 
liquid around the glass transition temperature. This technique per- 
mits us to explore the potential energy surface and barrier distribu- 
tions without suffering the exponential slowing down at low tempera- 
ture that affects molecular dynamics simulations. In agreement with 
previous studies we observe a distinct change in behaviour around 
T = 0.45, close to the dynamical transition temperature of mode 
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coupling theory (MCT). Below this temperature the number of dif- 
ferent local minima visited by the system for the same number of 
KMC steps decreases by more than an order of magnitude. The mean 
number of atoms involved in each jump between local minima and 
the average distance they move also decreases significantly, and new 
features appear in the partial structure factor. Above T ~ 0.45 the 
probability distribution for the magnitude of the atomic displacement 
per KMC step exhibits an exponential decay, which is only weakly 
temperature dependent. 
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1 Introduction 



A complete microscopic theory to explain the complex phenomenology of 
supercooled liquids and glasses is still elusive, and is the subject of a con- 
siderable research effort [0]. Mode-coupling theory, based upon a continuum 
approach to the wave vector dependence of correlation functions, has proved 
quite successful for supercooled liquids at higher temperatures, before acti- 
vated dynamics must be accounted for explicitly pj-Q. For finite systems, 
such as clusters, it has recently been possible to establish detailed connections 
between structure, thermodynamics, dynamics and the underlying potential 
energy surface (PES) 0. In the context of glasses, it was Goldstein who 
first noted the crucial role of the PES ||^, and Stillinger and Weber who first 
analysed computationally the local minima ('inherent structures') of various 
bulk models ^j. By analysing transitions between inherent structures they 
were able to identify the slow barrier crossings and localised rearrangements 
at low temperatures anticipated by Goldstein 0. 

Computer simulation has played a significant role in the testing and 



development of models for supercooled liquids and glasses |0. Interest 
in more direct connections to the PES has recently increased, with Sas- 



try, Debenedetti and Stillinger |]TT| characterising 'landscape- influenced' and 
'landscape-dominated' regimes for a binary Lennard- Jones system, in agree- 
ment with the instantaneous normal modes picture of Donati, Sciortino and 
Tartaglia [0. A number of other studies have also begun to characterise 
features of the underlying PES in more detail [p!B|-pD|. However, it has not 
yet been possible to account for the properties of supercooled liquids and 
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glasses at the level of detail achieved for finite systems 0. 

One obvious problem with the simulation of low temperature super- 
cooled liquids by standard molecular dynamics (MD) techniques is that the 
exploration of the PES slows down exponentially with decreasing tempera- 
ture. In the present contribution we show how a kinetic Monte Carlo (KMC) 
approach |^,|22| can provide useful information on the dynamics of a binary 
Lennard- Jones system around the glass transition temperature. This tech- 
nique combines accurate calculations of elementary transition states and the 
associated pathways with approximate dynamics, and hence it can overcome 
low temperature trapping and provide access to much longer simulation time 
scales. Hence it provides a view complementary to mode-coupling theory and 
standard MD simulations, since the rate theory employed is most accurate at 
low temperatures, and loses its validity at high temperatures. The paper is 
organised as follows. In section II we introduce the model and discuss the rel- 
evant points of the KMC method used. In sections III and IV we present the 
details of the KMC approach that we have used, and the numerical results, 
and in the last section we finish with the conclusions. 



2 The Model 

The system we have studied is a binary mixture of = 60 atoms in a cubic 
box with periodic boundary conditions, where 80% of the atoms are of type 
A and 20% are of type B. They interact via a Lennard- Jones potential of 
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the form 



Tap 
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Tap 



where r^p is the distance between particles and a, (3 G {A,B}. The param- 
eters we used were eab = 1-5, ess = 0.5, aAB = 0.8 and asB = 0.88, and 
all quantities will be measured in reduced units: energy in eAA-, length in 
aAAi temperature in tAA/kB {kB is the Boltzmann constant) and time in 
(a^^m/e^^)-'^/^, with m the mass of both types of atom. The density was 
fixed at 1.2. We employed the Stoddard-Ford shifted, truncated version of 
the above potential [^], with a cutoff at = 1.842, along with the minimum 
image convention |24]. 

Binary systems such as this have proved very popular amongst the 
glass simulation community, since they do not crystallise on the molecular 
dynamics (MD) time scale [11, 12 ,|l^, 18, 19, 25-37]. Stillinger and Weber were 
probably the first to consider an 80:20 mixture in their simulation of Ni8oP2o 



|38|. The present parameterisation was introduced by Kob and Andersen |25 
when they modified the potential of Ernst et al. because they found 
that it crystallised at low temperatures. The free energy barrier separating 
supercooled liquid from crystal is sufficiently high for the present parameters 
that crystallisation has not been reported in conventional MD studies |JD|. 

Following Angell 1^-^, glass-forming liquids can be classified as strong 
or fragile. Fragile systems exhibit non-Arrhenius temperature dependence of 
transport properties such as the diffusion constant, usually accompanied by 
a significant heat capacity peak at the glass transition. In contrast, strong 
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systems exhibit Arrhenius dynamics and small or negligible changes in ther- 
modynamic properties at the glass transition. Previous work has shown that 
the present system exhibits a significant degree of non- Arrhenius behaviour 
at low temperatures p5|-[27[]. 

Although the supercell size of only 60 atoms is rather small it has been 
shown in previous work to retain many of the most important features of the 



supercooled state [T^,|T3,|19|,^. Our aim in the present work is to analyse 
KMC trajectories consisting of a fixed number of elementary rearrangements 
between local minima as a function of temperature. In this preliminary study 
we report some basic structural and thermodynamic results, for comparison 
with previous work, and focus in particular on the changes that occur in the 
rearrangements sampled at lower temperatures. The energies of stationary 
points will be reported per atom, but energy barriers will be reported per 
supercell, since they scale intensively, not extensively, with the supercell size. 



3 The Kinetic Monte Carlo Approach 

Within the landscape picture we are concerned with the set of all particle 
positions, rj with i = 1,2,... ,N, which we can view as a 3A^-dimensional 
vector R, and its hopping dynamics on the PES, V^(R). For our purposes 
the most important points on the PES are the minima and the transition 
states that connect them. Here we follow Murrell and Laidler and define a 
true transition state as a saddle point with Hessian index one |^5[, i.e. the 
Hessian has only one negative eigenvalue whose eigenvector corresponds to 
the reaction coordinate and connects two minima. True transition states are 
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expected to dominate the dynamics in the low temperature regime because of 
the Murrell-Laidler theorem: if two minima are connected by a higher index 
saddle then there must be a lower energy path mediated by one or more true 
transition states |^ . 



The hopping dynamics between minima were simulated using a KMC 
approach |2^. We started from an initial random configuration and 
searched for a local minimum using Nocedal's limited memory LBFGS rou- 
tine [0. From the current local minimum we searched for transition states 



that connect it to other local minima. These searches were performed using 



a hybrid BFGS/eigenvector- following algorithm |^ where diagonalisation 
of the Hessian matrix is avoided. In fact, the Hessian is not required at all, 
but for the present empirical potential it proved more efficient to calculate 
the uphill search direction from the Hessian than to use a variational ap- 
proach To find more than one transition state for each minimum we 
used the method introduced by Middleton and Wales, where starting points 
for transition state searches are proposed by considering hard-sphere-type 



moves pO|. We employed a maximum of twenty transition states searches 
from each minimum. The minima connected to each transition state were 
then obtained from approximate steepest-descent paths commencing paral- 
lel and antiparallel to the transition vector using LBFGS energy minimisa- 
tion . All stationary points were finally converged to a root-mean-square 
force of less than 10~^ using one or two eigenvector-following steps, employ- 
ing full diagonalisation of the analytic Hessian matrix, to ensure that the 
stationary points have the correct Hessian index. 

The rate constants, used to calculate the probability per unit time of 



jumping to a new minimum, were calculated using Rice-Ramsperger-Kassel- 
Marcus (RRKM) theory within the harmonic approximation P5|,^: 

3N-3 

a n < 

^1^' ~ 3iV-4 ^ ' 

n ^ 

a=l 

where gi and gij are the order of the symmetry group of minimum i and tran- 
sition state ij (which are both almost always unity for the present system); 

and z/^^ are the frequencies of minimum i and transition state ij\ AE is 
the potential energy difference between minimum i and the transition state, 
and T is the temperature. 

If all the rate constants out of minimum i are known then we can 
evaluate the jump probability between two minima i and j as 

k 

where the sum is over all transition states. In the KMC method we randomly 
select one of the connected minima according to the probabilities defined in 
Eq. (^. Finally, once the system has moved to a new minimum the process 
is repeated. 

As in previous practical implementations of the KMC approach [Q, a 
further approximation is implicit because the sampling of transition states 
connected to minimum i is likely to be incomplete. However, this error does 
not affect the relative probabilities among the pathways located, but only 
the estimate of the corresponding time scale, which is approximate in any 
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case because of the model rate constants employed. Furthermore, we expect 
the omitted transition states to correspond to higher barrier processes with 
smaller rate constants on the basis of statistics obtained from more extensive 
surveys of stationary points . The application of KMC also carries with 
it the assumption that the jumps between minima are uncorrelated, so that 
the system equilibrates for long enough in each minimum to lose memory of 
the previous trajectory. We expect this assumption to be well satisfied in the 
low temperature regime below about T = 0.5, which is our main concern in 
the present work. As we noted above, the low temperature regime in which 
the KMC approach is most accurate should be complementary to approaches 
such as mode-coupling theory and standard simulation methods. 

The KMC approach differs in several ways from the activation-relaxation 



technique (ART) of Barkema and Mousseau et al. ||5T|-p6| . We calculate the 
transition states and pathways to very high precision, and locate more than 
one transition state per minimum in order to implement the KMC procedure. 
In previous ART studies the objective has generally been to achieve efficient 
structural relaxation, not to simulate true dynamics, and a Metropolis crite- 
rion has been used for jumps between minima using only the potential energy 
difference, without regard for effects such as vibrational entropy. 



4 Results 

We have performed a series of numerical KMC simulations to compute differ- 
ent properties of the liquid as we vary the temperature. We used a random 
starting configuration in each case, and hence it is important to discard the 
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initial part of the trajectory when the system is still equilibrating. In Figure 
|I] we plot the running average energy per atom of the local minima at KMC 
step n, defined as 

n 

(^)n=-E^- (4) 

where Vj is the energy per atom of the local minimum at KMC step j. In 
the course of a simulation, and before averaging other quantities, we have to 
be sure that at least local equilibration has been achieved. We assume that 
the system has equilibrated sufficiently when {V)^ reaches a stable value, 
which occurs after about 5000 KMC steps for all temperatures (Figure p. 
However, for the lowest temperature runs it is obvious that complete equi- 
libration has not been achieved, since the crystal should be the lowest free 



energy minimum [^. On the other hand, at low temperature the system 
samples fewer minima for KMC runs consisting of equal numbers of steps, 
and consequently the local equilibration actually appears to be faster. 

Our mean energies, averaged over the last 5000 KMC steps, are simi- 
lar to those obtained in reference using MD simulation (and a slightly 



different cut-off), except for our value at T = 1, which appears to be signif- 
icantly too low. Similar behaviour has been observed before in comparison 
of ART and MD results for silica, where the local minima obtained by ART 
sampling at the highest temperature were rather lower than those obtained 
from MD. For the KMC case the failure at high temperature is likely due 
to the breakdown in the Markov assumption for the minimum to minimum 
model dynamics. The system does not stay long enough in each minimum to 
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establish equilibrium properties and lose its memory of previous jumps. The 
transition rate between high energy minima is therefore probably underesti- 
mated, and a systematic bias towards lower energy minima results. 

The regions of the PES sampled in the present work are all high in 
potential energy compared to the lowest minimum that we found for this 
system, which is shown in Figure In this crystalline structure, based on 
a distorted face-centred-cubic lattice, the B atoms are arranged in a zig-zag 
pattern This minimum (energy —4.71 per atom) was found using a new 



stochastic global optimisation method, which combines the basin-hopping 
approach with KMC-type steps between local minima; further details 

will be published elsewhere. 

A qualitative change in behaviour at low temperature is clearly seen 
in Figure |^, where we plot the number of different local minima (including 
permutational isomers) visited in the last 5000 KMC steps of each run. The 
number of different local minima visited declines steadily with temperature 
until around T = 0.45, when it drops by more than an order of magnitude. 

In order to elucidate whether the system is a supercooled liquid or a 
glass, we calculated the partial structure factor, which encodes information 
about the typical distances between particles and their spatial distribution 
in a given minimum: 

\i=l j=l I 

where the average is over different local minima and fifty different orientations 
of q uniformly distributed over the surface of a sphere. We have calculated 
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Sajsio) for the A atoms over a range of temperatures (Figure averaging 
over every 25th local minimum from the last 5000 minima of each KMC 
simulation. There appears to be little dependence of Sapi^q) on temperature 
above T = 0.45, but new peaks appear at lower temperature and Sap{q) 
attenuates more slowly as a function of q. These results are in good agreement 
with previous work [p^^. 

We have also calculated the variance of the potential energy of the last 
5000 local minima for each run, as a first approximation to the configurational 
component of the heat capacity. Figure ^ shows that this variance exhibits 
a peak around T = 0.45. Figure |^ shows the mean number of A atoms that 
move more than a distance r, again averaged over the final 5000 local minima 
of each KMC run: 



, 10000 Na 

^^W = H000 5Z Y.®{K{3)-r^{j + l)\-r), (6) 

j=5001 a=l 

where is the Heaviside step function and is the position vector of atom 
a in minimum j. Above T = 0.45 Na{t) exhibits an exponential decay, while 
at lower temperatures the behaviour is more complicated. Above T = 0.45 
we can define a characteristic length in terms of the inverse of the exponential 
decay rate, which seems to be only weakly dependent on T. 

Below T = 0.45 the absence of a clear exponential decay may be due 
to heterogeneities in the dynamics. The atomic displacements tend to be 
shorter and the number of atoms involved is also reduced for both A and B 
atoms. These results are in good agreement with Mousseau's ART results 



for a much larger binary Lennard- Jones system |5^ and with the experimen- 



12 



tal observations of Tang et al. for some metallic glasses. Non-ergodic 
behaviour is again evident from the irregular behaviour of the plots at low 
temperature. Only local equilibration is achieved at such temperatures for 
KMC runs of the present length, and so the results depend sensitively upon 
the limited region that the system then samples. Nevertheless, the sampling 
should be significantly better than would be achieved by conventional MD. 

Figure |^ shows probability distributions for the barriers overcome in the 
last 5000 KMC steps as a function of temperature. The distributions were 
constructed in the same way as reference [^, and provide an interesting com- 
parison, since the canonical KMC sampling is different from the scheme used 
by Middleton and Wales. The distributions peak at low energy for the higher 



temperature runs, in agreement with the previous samples ||20| . However, for 
the two lowest temperatures there is a distinct shift in the maximum of the 
distribution to higher energy. The reason for this discrepancy has not yet 
been elucidated, but the new results suggest that the facile rearrangements 
most often encountered at higher temperature are no longer available at the 
lowest temperatures. However, a subsidiary maximum in the distribution 
is clearly visible around a barrier height of 0.8 at T = 0.4 and T = 0.425. 
We have previously found that such peaks correspond to diffusive-like pro- 
cesses in a number of model glasses where at least one atom changes 
its nearest-neighbour coordination shell. For KMC steps where the distance 
between the local minima is greater than unity, the forward barrier is always 
found to exceed 1.5 reduced units. The main features of the present results 
are therefore in reasonable agreement with previous work. 
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5 Conclusions 



We have performed numerical simulations of a supercooled binary Lennard- 
Jones liquid using a kinetic Monte Carlo (KMC) approach, which enables us 
to perform more extensive sampling at low temperature than the standard 
molecular dynamics approach. Most previous simulation studies have re- 
ported some sort of qualitative change in behaviour as the system approaches 
the critical temperature of mode-coupling theory Tc ~ 0.435 p5| , p6| . 

Following Jonsson and Anderson |^T|, Sastry, Debenedetti and Stillinger |11[] 
found evidence that around T = 0.45, above the glass transition, the height of 
the barriers separating local minima increases abruptly, and that the system 
shows evidence of activated dynamics at about the same temperature. They 
referred to the onset of activated dynamics as the 'landscape-dominated' 
regime, in agreement with the results of Donati, Sciortino and Tartaglia, 
based on an instantaneous normal modes picture |]12|. Sastry, Debenedetti 
and Stillinger also noted the onset of non-exponential relaxation below about 
T = 1.0, and referred to this region as 'landscape-influenced'. Local minima 
obtained by quenching from configurations generated at T = 0.5 were also 
found to escape to different local minima much more easily than local min- 
ima obtained from configurations generated at T = 0.4. The latter result is 
in good agreement with the change in the number of different local minima 
sampled in our calculations, suggesting that we have indeed sampled suffi- 
cient transition states per minimum to extract meaningful dynamics. This 
aspect of the KMC calculations will be investigated more systematically in 
future work; we intend to check that the results are converged with respect to 
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the number of transition state searches and are not affected by any samphng 
bias from the search algorithm. 

For a BLJ mixture of 50:50 composition Schr0der et al. concluded that 
the liquid dynamics could be separated into transitions between minima and 
vibrational motion, as originally suggested by Goldstein 0], below a temper- 
ature around Tc . They also noted correlated motions of groups of atoms 
in the same temperature regime, in agreement with results for an 80:20 BLJ 



mixture p8| , which suggest that the number of atoms in the correlated groups 
grows with decreasing temperature. This correlated motion is not inconsis- 
tent with our finding that the number of atoms moving in an elementary step 
decreases with temperature, because it involves sequences of such steps. 

Most recently, two groups have reported that the typical Hessian index 
of stationary points sampled by BLJ systems extrapolates to zero, again 
around T^. ||TB|,0. Our results show that transition states are still accessible 



below Tc, but support the general view that the PES sampled by the BLJ 
system changes in character somewhere around this point. 

In the high temperature regime above about T = 0.45 the system sam- 
ples the basins of attraction of many different local minima, and the distribu- 
tion of atomic displacements per KMC step exhibits exponential decay with 
a relatively weak temperature dependence. Hence we can define a character- 
istic length scale based upon the inverse of the exponential decay rate. In the 
low temperature regime the system samples far fewer local minima for the 
same number of KMC steps, and the atomic dynamics cannot be described 
by a single length scale. We interpret this phenomenon in terms of a multi- 
funnel energy landscape [^, where the system is associated with lower lying 
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minima at lower temperature leading to a higher effective activation energy 
for diffusion, and non-Arrhenius dynamics p2|-p5|]. 

It is noteworthy that this popular binary Lennard- Jones system does 
in fact possess a crystalline phase well separated in energy from the minima 
sampled by the supercooled liquid [0. In fact, crystallisation was previously 
observed by Kob and Anderson p6[ for a system of the same composition but 
with different e parameters. The present approach can also be used to study 
dynamical properties such as self-diffusion or the incoherent intermediate 
scattering function. This work is now in progress. 
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Figure 1: Running average of the potential energy per atom, (V)^, as a 
function of the number of KMC steps, n, at temperatures T — 1, 0.6, 0.5, 
0.45 and 0.425. We indicate with an arrow the energy of the lowest minimum 
that we have found. 
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Figure 2: Lowest minimum found for the binary Lennard- Jones mixture. 
The 48 A atoms and 12 B atoms are shaded gray and black, respectively. 
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Figure 3: Number of different local minima visited in the last 5000 KMC 
steps as a function of the temperature. 
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Figure 4: Partial structure factor for A atoms as a function of temperature. 
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Figure 5: Mean square fluctuation of tlie potential energy divided by as 
a function of T. 



26 



100 




0.0 0.1 0.2 0.3 



r 



Figure 6: Linear- log plot of the mean number of A atoms moving through 
a distance r, NA{r), for KMC runs at five temperatures. 
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Figure 7: Normalised barrier distributions for accepted KMC moves as a 
function of temperature. Only the barriers corresponding to forward jumps 
(not backward) are included. 
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